* ---------------------------------------------
* ' test individual vs contextual vs minority-majority index
* ' across different levels
* ---------------------------------------------
global home_dir "/N/project/suicide_study/pnas_replication"

args model model_version 

local geo_type "county"
*local model "mi"
*local model_version 6


global figure_file "${home_dir}/results/log"

cd ${home_dir}

local output_file "${figure_file}/`model'_`model_version'.smcl"

log using `output_file', replace 

if ("`model'" == "logit"){
	use "${home_dir}/data/processed/suicide_reg_v1_raw.dta", clear
}

if ("`model'" == "mi") {
	use "${home_dir}/data/processed/suicide_reg_v1_imputed_M10.dta", clear 	
}

* for now, we use the following simple survey weights 
if ("`model'" == "mi") {
	mi svyset `geo_type' [pw=ObsWgt0] 
}
else {
	svyset `geo_type' [pw=ObsWgt0] 	
}

* margins for each category
program margin_interact 
	args X Y k model
	sum `X', d 
	local gap = (`r(max)' - `r(min)') / `k' 
	if ("`model'" == "logit") {
		margin `Y', at(`X' = (`r(min)' (`gap') `r(max)')) predict(pr)
	} 
	else if ("`model'" == "mi") {
		mimrgns `Y', at(`X' = (`r(min)' (`gap') `r(max)')) predict(pr)
	}	
end 

program mchange_mi
	args X k model
	if ("`model'" == "logit") {

		if ("`k'" == "continuous") {
			sum `X' if e(sample), d 
			margin, at(`X' = (`r(min)' `r(max)')) post predict(pr)
			mlincom  2 - 1, decimal(7) stat(all)	
		} 
		else if ("`k'" == "binary") {
			sum `X' if e(sample), d 
			margin, at(`X' = (`r(min)' `r(max)')) post predict(pr)
			mlincom  2 - 1, decimal(7) stat(all)
		} 
		else if ("`k'" == "categorical") {
			margin `X' if e(sample)==1 , at() pwcompare predict(pr)	
		}
	}

	else if ("`model'" == "mi") {

		if ("`k'" == "continuous") {
			sum `X' , d 
			mimrgns, at(`X' = (`r(min)' `r(max)')) post predict(pr)
			mlincom  2 - 1, decimal(7) stat(all)
		} 
		else if ("`k'" == "binary") {
			sum `X' , d 
			mimrgns, at(`X' = (`r(min)' `r(max)')) post predict(pr)
			mlincom  2 - 1, decimal(7) stat(all)
		} 
		else if ("`k'" == "categorical") {
			mimrgns `X' , at()   pwcompare predict(pr)	
		}
	}
end

* create some dummy codings
tab Race5, gen(race5_nh)
	rename race5_nh1 White_nh 
	rename race5_nh2 Black_nh 
	rename race5_nh3 AIAN_nh 
	rename race5_nh4 AsPI_nh 
	rename race5_nh5 Hispanic

tab MarStat5, gen(ms)
	rename ms1 Marrd5
	rename ms2 Widow5
	rename ms3 Divor5
	rename ms4 Separ5
	rename ms5 NvMar5

tab AgeGrp4, gen(ag) 
	rename ag1 Age_15_24
	rename ag2 Age_25_44
	rename ag3 Age_45_64
	rename ag4 Age_65_Up

destring St, replace 

* set-up equations
local religion Rat_GC_ProE Rat_GC_ProM Rat_GC_ProB Rat_GC_Cath Rat_GC_Jew Rat_GC_Oth
local contextual_control Rat_Poverty Rat_Mig_Cum Pop_Den

local religion Rat_GC_ProE Rat_GC_ProM Rat_GC_ProB Rat_GC_Jew Rat_GC_Oth
local contextual_control Rat_Poverty Rat_Mig_Cum Pop_Den

local demographics_raw i.Female c.RAT_Female i.AgeGrp4 c.RAT_AgeGrp4_2 c.RAT_AgeGrp4_3 c.RAT_AgeGrp4_4 i.Race5 c.RAT_Race5_2 c.RAT_Race5_3 c.RAT_Race5_4 c.RAT_Race5_5 i.BornUSA c.RAT_BornUSA i.MarStat5 c.RAT_MarStat5_2 c.RAT_MarStat5_3 c.RAT_MarStat5_4 c.RAT_MarStat5_5
local demographics_individual i.Female i.AgeGrp4 i.Race5 i.BornUSA i.MarStat5 
local demographics_county c.RAT_Female c.RAT_AgeGrp4_2 c.RAT_AgeGrp4_3 c.RAT_AgeGrp4_4  c.RAT_Race5_2 c.RAT_Race5_3 c.RAT_Race5_4 c.RAT_Race5_5 c.RAT_BornUSA c.RAT_MarStat5_2 c.RAT_MarStat5_3 c.RAT_MarStat5_4 c.RAT_MarStat5_5
local demographics_same i.Female c.std_same_prop_Sex i.AgeGrp4 c.std_same_prop_AgeGrp4 i.Race5 c.std_same_prop_Race5 i.BornUSA c.std_same_prop_BornUSA i.MarStat5 c.std_same_prop_MarStat5 
local demographics_inter i.Female##c.std_same_prop_Sex i.AgeGrp4##c.std_same_prop_AgeGrp4 i.Race5##c.std_same_prop_Race5 i.BornUSA##c.std_same_prop_BornUSA i.MarStat5##c.std_same_prop_MarStat5


if (`model_version' == 1){
	local model_eq i.Year `demographics_raw' `contextual_control' `religion' 
	local margin_demographics Female RAT_Female AgeGrp4 RAT_AgeGrp4_2 RAT_AgeGrp4_3 RAT_AgeGrp4_4 Race5 RAT_Race5_2 RAT_Race5_3 RAT_Race5_4 RAT_Race5_5 BornUSA RAT_BornUSA MarStat5 RAT_MarStat5_2 RAT_MarStat5_3 RAT_MarStat5_4 RAT_MarStat5_5 
}
if (`model_version' == 2){
	local model_eq i.Year `demographics_raw' `contextual_control' `religion' i.UnEmpl c.RAT_UnEmpl i.PhysProb c.RAT_PhysProb
	local margin_demographics Female RAT_Female AgeGrp4 RAT_AgeGrp4_2 RAT_AgeGrp4_3 RAT_AgeGrp4_4 Race5 RAT_Race5_2 RAT_Race5_3 RAT_Race5_4 RAT_Race5_5 BornUSA RAT_BornUSA MarStat5 RAT_MarStat5_2 RAT_MarStat5_3 RAT_MarStat5_4 RAT_MarStat5_5 UnEmpl RAT_UnEmpl PhysProb RAT_PhysProb
}
if (`model_version' == 3){
	local model_eq i.Year `demographics_same' `contextual_control' `religion' 
	local margin_demographics Female std_same_prop_Sex AgeGrp4 std_same_prop_AgeGrp4 Race5 std_same_prop_Race5 BornUSA std_same_prop_BornUSA MarStat5 std_same_prop_MarStat5 
}
if (`model_version' == 4){
	local model_eq i.Year `demographics_same' `contextual_control' `religion' i.UnEmpl c.std_same_prop_UnEmpl i.PhysProb c.std_same_prop_PhysProb
	local margin_demographics Female std_same_prop_Sex AgeGrp4 std_same_prop_AgeGrp4 Race5 std_same_prop_Race5 BornUSA std_same_prop_BornUSA MarStat5 std_same_prop_MarStat5 UnEmpl std_same_prop_UnEmpl PhysProb std_same_prop_PhysProb
}
if (`model_version' == 5){
	local model_eq i.Year `demographics_inter' `contextual_control' `religion' 
}
if (`model_version' == 6){
	local model_eq i.Year `demographics_inter' `contextual_control' `religion' i.UnEmpl##c.std_same_prop_UnEmpl i.PhysProb##c.std_same_prop_PhysProb
}	

if (`model_version' == 7){
	local model_eq i.Year `demographics_individual' `contextual_control' `religion' 
	local margin_demographics Female AgeGrp4 Race5 BornUSA MarStat5 
}
if (`model_version' == 8){
	local model_eq i.Year `demographics_county' `contextual_control' `religion' 
	local margin_demographics RAT_Female RAT_AgeGrp4_2 RAT_AgeGrp4_3 RAT_AgeGrp4_4 RAT_Race5_2 RAT_Race5_3 RAT_Race5_4 RAT_Race5_5 RAT_BornUSA RAT_MarStat5_2 RAT_MarStat5_3 RAT_MarStat5_4 RAT_MarStat5_5 
}

if (`model_version' == 9){
	local model_eq i.Year `demographics_individual' `contextual_control' `religion'  i.UnEmpl i.PhysProb
	local margin_demographics Female AgeGrp4 Race5 BornUSA MarStat5 UnEmpl PhysProb
}
if (`model_version' == 10){
	local model_eq i.Year `demographics_county' `contextual_control' `religion'  c.RAT_UnEmpl c.RAT_PhysProb
	local margin_demographics RAT_Female RAT_AgeGrp4_2 RAT_AgeGrp4_3 RAT_AgeGrp4_4 RAT_Race5_2 RAT_Race5_3 RAT_Race5_4 RAT_Race5_5 RAT_BornUSA RAT_MarStat5_2 RAT_MarStat5_3 RAT_MarStat5_4 RAT_MarStat5_5 RAT_UnEmpl RAT_PhysProb
}

* test how long it would take.
* mi estimate: svy: mean Suic 
* local demographics i.Female c.RAT_Female i.AgeGrp4 c.RAT_AgeGrp4_2 c.RAT_AgeGrp4_3 c.RAT_AgeGrp4_4 i.Race5 c.RAT_Race5_2 c.RAT_Race5_3 c.RAT_Race5_4 c.RAT_Race5_5 i.BornUSA c.RAT_BornUSA i.MarStat5 c.RAT_MarStat5_2 c.RAT_MarStat5_3 c.RAT_MarStat5_4 c.RAT_MarStat5_5
* mi estimate: svy: logit Suic i.St `demographics' UnEmpl RAT_UnEmpl PhysProb RAT_PhysProb

* main effects : margins
if ("`model'" == "mi"){

	mi estimate: svy: logit Suic i.St `model_eq', or 
	estimates store m1 

	if (`model_version' <= 4 | `model_version' >= 9){
		if (`model_version' != 10){
			estimates restore m1
			mchange_mi Female "binary" "mi"
			estimates restore m1
			mchange_mi AgeGrp4 "categorical" "mi"
			estimates restore m1
			mchange_mi Race5 "categorical" "mi"
			estimates restore m1
			mchange_mi BornUSA "binary" "mi"
			estimates restore m1
			mchange_mi MarStat5 "categorical" "mi"
		}
		
		if (`model_version' == 1 | `model_version' == 2 | `model_version' == 10) {
			estimates restore m1
			mchange_mi RAT_Female "continuous" "mi"
			estimates restore m1
			mchange_mi RAT_AgeGrp4_2 "continuous" "mi"
			estimates restore m1
			mchange_mi RAT_AgeGrp4_3 "continuous" "mi"
			estimates restore m1
			mchange_mi RAT_AgeGrp4_4 "continuous" "mi"
			estimates restore m1
			mchange_mi RAT_Race5_2 "continuous" "mi"
			estimates restore m1
			mchange_mi RAT_Race5_3 "continuous" "mi"
			estimates restore m1
			mchange_mi RAT_Race5_4 "continuous" "mi"
			estimates restore m1
			mchange_mi RAT_Race5_5 "continuous" "mi"
			estimates restore m1
			mchange_mi RAT_BornUSA "continuous" "mi"
			estimates restore m1
			mchange_mi RAT_MarStat5_2 "continuous" "mi"
			estimates restore m1
			mchange_mi RAT_MarStat5_3 "continuous" "mi"
			estimates restore m1
			mchange_mi RAT_MarStat5_4 "continuous" "mi"
			estimates restore m1
			mchange_mi RAT_MarStat5_5 "continuous" "mi"
		}
		if (`model_version' == 3 | `model_version' == 4) {
			estimates restore m1
			mchange_mi std_same_prop_Sex "continuous" "mi"
			estimates restore m1
			mchange_mi std_same_prop_AgeGrp4 "continuous" "mi"
			estimates restore m1
			mchange_mi std_same_prop_Race5 "continuous" "mi"
			estimates restore m1
			mchange_mi std_same_prop_BornUSA "continuous" "mi"
			estimates restore m1
			mchange_mi std_same_prop_MarStat5 "continuous" "mi"
		}
	
		if (`model_version' == 2 | `model_version' == 4 | `model_version' == 9 ){
			estimates restore m1
			mchange_mi UnEmpl "categorical" "mi"
			estimates restore m1
			mchange_mi PhysProb "categorical" "mi"
		}
	
		if (`model_version' == 2 | `model_version' == 10){
			estimates restore m1
			mchange_mi RAT_UnEmpl "continuous" "mi"
			estimates restore m1
			mchange_mi RAT_PhysProb "continuous" "mi"
		}
		if (`model_version' == 4){
			estimates restore m1
			mchange_mi std_same_prop_UnEmpl "continuous" "mi"
			estimates restore m1
			mchange_mi std_same_prop_PhysProb "continuous" "mi"
		}		
	}
}

if ("`model'" == "logit"){
	svy: logit Suic i.St `model_eq', or 
	estimates store m1 

	if (`model_version' <= 4 | `model_version' == 7 | `model_version' == 8){
		mchange `margin_demographics', amount(all) delta(100) statistics(all) decimals(7)
	}

}

if (`model_version' == 5 | `model_version' == 6) {
	estimates restore m1
	margin_interact std_same_prop_Sex Female 5 `model'
	estimates restore m1
	margin_interact std_same_prop_AgeGrp4 AgeGrp4 5 `model'
	estimates restore m1
	margin_interact std_same_prop_Race5 Race5 5 `model'
	estimates restore m1
	margin_interact std_same_prop_BornUSA BornUSA 5 `model'
	estimates restore m1
	margin_interact std_same_prop_MarStat5 MarStat5 5 `model'

	if (`model_version' == 6) {
		estimates restore m1
		margin_interact std_same_prop_UnEmpl UnEmpl 5 `model'
		estimates restore m1
		margin_interact std_same_prop_PhysProb PhysProb 5 `model'
	}

}


log close 
